#ifndef __tet_hh__
#define __tet_hh__

class Tet
{
public:
  Tet(const double t[3*4]) : vertices(t) { compute_volume() ; }
  virtual ~Tet() {}
  virtual double *evaluate_basis ( int n_points, const double points[] ) = 0 ;
  virtual int test() = 0 ;
  
  virtual double* reference_to_physical_coords ( int n_points, double ref[] ) ;
  virtual double* physical_to_reference_coords ( int n_points, double phy[] ) ;
  
  double* barycentric_coords ( double p[3] ) ;
  virtual double element_volume () ;
  
protected:
  const double* vertices ;
  double volume;

private:
  void compute_volume() ;

};

inline
void Tet::compute_volume()
{
//
//           | x1 x2 x3 x4 |
//  Volume = | y1 y2 y3 y4 |
//           | z1 z2 z3 z4 |
//           |  1  1  1  1 |
//
// Alternately, we can directly compute the determinant:
//
//  double a[4*4];
//  int i, j;
//
//  for ( i = 0; i < 3; i++ )
//    for ( j = 0; j < 4; j++ )
//      a[i+j*4] = this->vertices[i+j*3];
//
//  i = 3;
//  for ( j = 0; j < 4; j++ )
//    a[i+j*4] = 1.0;
//
//  volume = fabs ( r8mat_det_4d ( a ) ) / 6.0;
//
  volume =
      vertices[0+0*3] * (
        vertices[1+1*3] * ( vertices[2+2*3] - vertices[2+3*3] )
      - vertices[1+2*3] * ( vertices[2+1*3] - vertices[2+3*3] )
      + vertices[1+3*3] * ( vertices[2+1*3] - vertices[2+2*3] ) )
    - vertices[0+1*3] * (
        vertices[1+0*3] * ( vertices[2+2*3] - vertices[2+3*3] )
      - vertices[1+2*3] * ( vertices[2+0*3] - vertices[2+3*3] )
      + vertices[1+3*3] * ( vertices[2+0*3] - vertices[2+2*3] ) )
    + vertices[0+2*3] * (
        vertices[1+0*3] * ( vertices[2+1*3] - vertices[2+3*3] )
      - vertices[1+1*3] * ( vertices[2+0*3] - vertices[2+3*3] )
      + vertices[1+3*3] * ( vertices[2+0*3] - vertices[2+1*3] ) )
    - vertices[0+3*3] * (
        vertices[1+0*3] * ( vertices[2+1*3] - vertices[2+2*3] )
      - vertices[1+1*3] * ( vertices[2+0*3] - vertices[2+2*3] )
      + vertices[1+2*3] * ( vertices[2+0*3] - vertices[2+1*3] ) );

  if ( volume <= 0.0 )
  {
    cerr << "\n";
    cerr << "BASIS_MN_TET4 - Fatal error!\n";
    cerr << "  Element has zero volume.\n";
    exit ( 1 );
  }
}

inline
double Tet::element_volume()
{
  return this->volume ;
}


class Tet4 : public Tet
{
public:
  Tet4(const double t[3*4]) : Tet(t) {} 
  virtual double *evaluate_basis ( int n_points, const double points[] ) ;
  virtual int test() ;
  
private:
};


class Tet10 : public Tet
{
public:
  Tet10(const double t[3*4]) : Tet(t) {} 
  virtual double *evaluate_basis ( int n_points, const double points[] ) ;
  virtual int test() ;
  
private:
};



#endif // #ifndef __tet_hh__
